function zero = sdp_solve_w(wbar)

global par_;

chi  = par_.chi;
menucost =  par_.menucost; 
std_eA = par_.std_eA;
alpha  = par_.alpha;
alphapos  = par_.alphapos;
phiR   = par_.phiR;
phiA   = par_.phiA;
std_m  = par_.std_m;
omega  = par_.omega;

% Unpack parameters
%mu_ = par_.mu_;
model = par_.model;
mfreq = par_.mfreq;
estim = par_.estim;
show = par_.show;

params_set;
matrices_set;

if lbar>1 || lbar<0 || rho<0 || rho>1 || menucost<0 || std_eA < 0
    zero = NaN; % tells solver to check another point
    disp('Limits reached on parameters')
    return;
end

Cbar = (wbar/chi)^(1/gam);  

PAYOFFMAT = Cbar*(PMAT.^(1-theta) - wbar*PMAT.^(-theta/(1-eta)));

v_iter;  

if exitflag==1
    
    Pdist=zeros(nump,1);    
    
    Pdist(ceil(nump/2),:) = 1;
   
    Pdist=Pdist/(sum(sum(Pdist)));   
   
    PdistDIFF=inf;           % reset difference of Pdist
    Piter=0;                 % reset counter
    while (PdistDIFF>sqrt(eps) && Piter<1000)  % iterate to convergence 
        Piter=Piter+1;                    % increment counter
        Pdist_sh = R*T2*Pdist;            % distribution after shocks
        PdistNew = (1-lambda).*Pdist_sh + P.*(ones(nump,nump)*(lambda.*Pdist_sh)); %P is logitprobMAT
        PdistDIFF=nump*max(max(abs(PdistNew-Pdist)));   % sup norm 
        Pdist=PdistNew;                                      % updating Pdist   
    end
    if PdistDIFF>sqrt(eps)
       disp('Pdist not converged!')
       sound(sin(1:300))
    end
    Pdist=Pdist/(sum(sum(Pdist)));   

    calcstats;

    par_.PD  = PD;
    par_.Cbar = Cbar;
    par_.freqpchanges  = freqpchanges;
    par_.MeanAbsPchange = MeanAbsPchange;
    par_.KurtPchanges = KurtPchanges;
        
    par_.lambda = lambda;
    par_.lambda_prime = lambda_prime;
    par_.Pdistans = Pdistans;
    par_.Pdist_sh = Pdist_sh;
        
    par_.Pchangegrid = Pchangegrid;
    par_.EPchange = EPchange;
    par_.StdPchanges = StdPchanges;
    par_.denspchanges = denspchanges;

    par_.Rss = Rss; 
    par_.Cbar = Cbar; 
    par_.phiPI = NaN; 
    par_.phiC = NaN; 

    par_.lbar = lbar; 
    %par_.ksi = ksi;              
    par_.noise = noise;

    par_.PMAT = PMAT; 
    par_.T2 = T2; 
    par_.Pgrid = Pgrid; 
    par_.pstep = pstep;
    par_.D = D;

    par_.Pdist = Pdist;
    par_.PD  = PD;
    par_.mbar = mbar;
    par_.Vf = Vf;

    par_.wbar = wbar;
    
    if show 
        printstats
        plotfigs
    else
        zero(1) = 1-sum(sum(PMAT.^(1-theta).*Pdist))^(1/(1-theta));
    end
    

    if  par_.dodyn
            switch par_.switch_sim
                case 'irf'
                    dynsim
                    plotirf
            end
      
    end
end
end
